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ABSTRACT 

A low-dimensional dynamical model for transitional buoyancy- driven flow in a differentially heated tall 
enclosure is presented. The full governing partial differential equations with the associated boundary con- 
ditions are solved by a spectral element method for a cavity of aspect ratio A ~ 20. Proper orthogonal 
decomposition is applied to the oscillatory solution at Prandtl number Pr — Pr Q = 0.71 and Grashof num- 
ber Gr — Gr 0 — 3.2 x 10 4 to construct empirical eigenfunctions. Using the four most energetic empirical 
eigenfunctions for the velocity and temperature as basis functions and applying Galerkin’s method, a reduced 
model consisting of eight nonlinear ordinary differential equations is obtained. Close to the “design” condi- 
tions (Pr 0 ,Gr 0 ), the low-order model (LOM) predictions are in excellent agreement with the predictions of 
the full model. In particular, the critical Grashof number at the onset of the first temporal flow instability 
(Hopf bifurcation) as well as the frequency and amplitude of oscillations at supercritical conditions are in 
excellent agreement with the predictions of the full model. Far from the “design” conditions, the LOM 
predicts the existence of multiple stable steady solutions at large values of Gr, and a unique stable steady 
solution at small values of Gr, and exhibits hysteretic behavior that is qualitatively similar to that observed 
in direct numerical simulations based on the full model. 

INTRODUCTION 

Thermally driven cavity flows provide a wealth of paradigms for the study of flow instabilities and transition 
to turbulence. A classical problem is concerned with the motion in a rectangular cavity when the vertical 
boundaries are maintained at fixed but distinct temperatures and the upper and lower horizontal surfaces 
are adiabatic. For a Newtonian fluid, subject to the Boussinesq approximation, the flow field is governed by 
three dimensionless parameters: the aspect ratio A (height/width), the Prandtl number Pr, and the Grashof 
number Gr (based on the cavity width) . At moderate Gr, the base flow corresponds to a single unicellular 
state. As Gr increases, both time-dependent and stationary bifurcations can occur, depending on the aspect 
ratio and Prandtl number. In the present study we investigate the possibility of developing low- dimensional 
dynamical models of oscillatory multicellular convection in a rectangular cavity of aspect ratio A = 20. 

The method of weighted residuals can be used to transform evolution partial differential equations (PDEs) to 
systems of ordinary differential equations (ODEs). Most frequently, trigonometric or orthogonal polynomials 
are used as basis functions (e.g., Gottlieb and Orszag, 1977). Other basis functions, e.g., spline functions 
(Liakopoulos, 1985, Liakopoulos and Hsu, 1984) have also been used successfully. In practice, the infinite 
dimensional problem is truncated to a finite dimensional space of dimension n. To ensure that the dynamic 
behavior described by the resulting finite dimensional system converges to that of the original infinite di- 
mensional problem, the required dimension n is typically high. A large reduction in the dimension of the 
resulting system of ODEs is accomplished by expanding the unknown functions in terms of basis functions 
that are constructed specifically for each flow system and reflect the behavior of the flow in the vicinity of 
some values of the controlling parameters. Proper Orthogonal Decomposition (POD) is a rigorous method- 
ology for obtaining a set of optimal basis functions (Berkooz et al., 1993). POD systematically identifies the 
most energetic eigenmodes that contain enough information for accurate description of the flow dynamics. 
Consequently, one is able to compress numerical or experimental data by retaining a small number of modes 
that capture most of the flow and temperature “energy”. 
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In this paper, we consider the buoyancy-driven flow in a laterally heated enclosure of aspect ratio A = 20. 
Direct numerical simulations are performed for Pr — 0.71 and 10 3 < Gr < 10 5 . Proper Orthogonal De- 
composition and Galerkin’s version of the method of weighted residuals are employed to derive low-order 
dynamical models of the flow in the transitional regime. 

FULL MODEL 

The flow domain consists of a two-dimensional rectangular cavity in which the upper and lower surfaces 
are adiabatic and the vertical walls are maintained at constant but different temperatures. Subject to the 
Boussinesq approximation, the continuity, momentum, and energy equations are written in dimensionless 
form: 

V ■ V = 0, (1) 


-jrr + (V • V)F + VP = 0j + -^=V 2 V 

ut y/Gr 

— + V • V0 = — 1 - 7 =v 2 q 
dt Pry/ Gr 


( 2 ) 


( 3 ) 


with boundary conditions V = 0 along the walls, 0(x = 0,y) = 1, 0(x = 1, y) — 0, and |~(x,±y) = 0. 

The dimensionless variables are defined as 


(x f y) = (** f y*)/Z, t= 

l u . 




p = and Q = 


P^c 


( 4 ) 


Here T\ and T 2 are the cold and hot wall temperatures, respectively, l is the cavity width, Pr — vja is 
the Prandtl number, Gr = is the Grashof number, v is the kinematic viscosity, a is the thermal 

diffusivity, /? is the thermal expansion coefficient, g is the acceleration of gravity, j is the unit vector in the 
direction opposite to gravity, and u c = y/(3gl(T2 — T\). Note that the characteristic velocity u c has been 
determined by balancing the inertia and buoyancy forces in the momentum equation and the pressure has 
been scaled by two times the dynamic pressure. These are appropriate scales for the high Grashof number 
thermal convection considered in the present study. 


For high Gr the IBVP (full model) has multiple stable solutions that can be computed by assigning 
appropriate initial conditions. For example, if we start with a three-cell time-independent solution at 
Gr = 4 X 10 4 and gradually decrease Gr, different solution branches can be reached. At Gr = 3.2 x 10 4 , 
a stable three-cell time-independent solution is found. This contrasts the four-cell time-dependent flow that 
is found when Gr is gradually increased. Reducing Gr further, we obtain a sequence of time-independent 
solutions that may belong to the same steady solution branch. At Gr = 10 4 a steady multicellular solution is 
computed. Note that the multiplicity of solutions is a property of the full model for large values of Gr only. 
At Gr = 5 X 10 3 , the solution to the IBVP computed by decreasing Gr is identical to the one computed 
for increasing Gr. We have found no evidence of multiple solutions as Gr was decreased further. Although 
an exhaustive search for all possible solutions is prohibitively expensive, we may conclude with a reasonable 
degree of confidence that the full model has a unique stable steady solution for Gr < 8 X 10 3 . 
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DERIVATION OF THE LOW-ORDER MODEL 


In the context of transitional thermal convection, we assume that, for some values ( Pr oy Gr 0 ), we have 
obtained spontaneously oscillatory flow and temperature fields. We refer to these values of the parameters 
as “design” parameters or “design” conditions. Furthermore, it is assumed that M snapshots of each field 
have been experimentally measured or computed based on the full model for Pr — Pr 0 andGr = Gr 0 . 
Following the procedure described in Liakopoulos and Gunes (1996), we construct the stationary empirical 
eigenfunctions. Expanding the unknown functions in terms of the eigenfunctions, applying Galerkin method, 
and making use of the orthonormality property of the empirical eigenfunctions, we obtain a system of non- 
linear ODEs for the temporal expansion coefficients 

^7“ ~ A/c + ~^=B k + Ck% a x + —=D kx a x + Ek t jCL x aj + R kx b x (5) 

at vGr VGr 

= Fk + Pr^/Gr ^ 1 * + Hkt a-x + p r l/o;Il + Jkij^j + K kx b t . 

EIGENVALUES AND EMPIRICAL EIGENFUNCTIONS 

Twenty snapshots (Af = 20) of the oscillatory solution for Gr — Gr 0 — 3.2 x 10 4 and Pr = Pr Q = 0.71 
provide the input data for the extraction of the empirical eigenfunctions. A list of the eight largest eigenvalues 
is given in Table I. The eigenvalues have been normalized so that their sum equals unity, and their cumulative 
contribution to the total flow and temperature fluctuation “energy” is given in the last column of Tables 
la and lb, respectively. The first four eigenmodes capture 99.83% of both flow and temperature fluctuation 
“energy”. However, the energy distribution among the four most energetic modes is slightly different for 
temperature than for velocity. The principal eigenvalue contributes 63.07% of the total fluctuation “energy” 
for temperature as opposed to 55.83% for velocity. Furthermore, for the temperature eigenvalues A 2 ~ A ; /2 
and A 4 ~ A 3 / 2 . 

Isotherms and streamlines for the four most energetic empirical eigenfunctions are shown in Figure 1. Note 
that although the direct numerical simulations and the proper orthogonal decomposition are performed 
in terms of primitive variables (u, u, P, 0), the velocity eigenfunctions are presented here in terms of the 
corresponding streamfunction. The velocity and temperature eigenfunctions are both centro-symmetric. 
The two most energetic eigenfunctions are localized in the central part of the cavity where velocity and 
temperature fluctuations are most vigorous. Higher order eigenfunctions capture features in the entire 
cavity, including regions close to the horizontal end walls. Repeated flow and temperature structures can 
be seen in the middle part of the cavity. These repeated structures form patterns that are characteristic of 
the early transition process of convection in laterally heated tall cavities. Preliminary results reveal that the 
most energetic eigenmodes exhibit the same spatial patterns in the middle part of a cavity of aspect ratio 
A = 40. 


RESULTS BASED ON THE LOW-ORDER MODEL 

In developing a reduced model, it is desired to keep the dimensionality of the system sufficiently low so that 
the methods of the dynamical systems theory can be effectively applied. On the other hand, enough modes 
should be retained so that the field variables are reconstructed accurately, most of the flow and tempera- 
ture fluctuation “energy” is captured, and the potentially important information hidden in the small scale 
features of higher modes is not lost. For the problem at hand, at least four modes for each field variable 
need to be retained in order to obtain oscillations of correct amplitude at “design” conditions ( Gr 0 , Pr 0 ). 
Retaining fewer than four modes leads to unrealistically large modal amplitudes. The results summarized 
below are obtained using eight ODEs (four from the momentum equation and four from the energy equa- 
tion). The critical points (fixed points, steady solutions) of the low-order model are found by setting the 
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right hand side of equation (5) to zero and solving the resulting nonlinear algebraic equations. The stability 
of a steady solution is then determined by the eigenvalues of the associated Jacobian matrix evaluated at 
the critical point under consideration. Note that the Jacobian matrix is real and non-symmetric. Several 
steady solution branches are found for 1< Gr < 2 X 10 5 using a Newton-Raphson method with random 
initial guesses for the solution components. For Gr < 2.63 X 10 4 , a unique branch of stable fixed points 
is found. We refer to this branch as the primary solution branch. It is denoted by A in Figure 2 where 
the norm (a^a* + 6A ) 1 ^ 2 of the fixed points is plotted against Gr. For Gr > 2.63 x 10 4 , multiple steady 
solutions are found. Since we are interested in solutions that are potentially relevant to the solutions of 
the original system of PDEs, we discard fixed points whose norm is several orders of magnitude larger than 
the norm of the primary branch. In Figure 2, we present only fixed points with a norm smaller than 3 in 
the range 1< Gr < 10 5 . Thick lines represent stable fixed points while thin lines correspond to unstable 
fixed points. Two stable steady solutions exist for 2.63 xlO 4 < Gr < 3.11 X 10 4 . At Gr = 3.11 X 10 4 , the 
fixed point on the primary branch undergoes a Hopf bifurcation that marks the onset of periodic oscilla- 
tions in time. For 3.11 xlO 4 < Gr < %.14 X 10 4 , we find one stable steady solution (branch B) while for 
4.14xl0 4 < Gr < 9.87 x 10 4 , two stable steady solutions exist (branches B and C). The fixed point along 
branch C becomes unstable at Gr = 9.87 X 10 4 . For Gr > 9.87 X 10 4 , we find only one stable steady solution 
that maintains its stability at least up to Gr — 2 X 10 5 . These stability results have been confirmed by solving 
equations (5) with suitable initial conditions using a fourth order Runge-Kutta ODE solver. 

A detailed quantitative comparison of the low-order model predictions with the full model is beyond 
the scope of the present work. However, we mention some encouraging observations both quantitative and 
qualitative in nature. The low-order model prediction of a Hopf bifurcation along the primary branch is 
in excellent agreement with the full model prediction. The low-order model predicts a Hopf bifurcation 
at Gr = Gr/f = 3.11 X 10 4 , while the full model predicts the onset of spontaneously oscillatory convection 
at Gr ~ 3.1 X 10 4 . The frequency of oscillations at the onset of the temporal instability is successfully 
predicted by the low-order model. Turning our attention to values of Gr far from the design conditions, 
we find qualitative agreement between the low-order model and the full model. For example, the low-order 
model predicts a stationary bifurcation at Gr < Gr/j in qualitative agreement with the full model which 
predicts that a stationary instability certainly precedes the onset of oscillatory convection for A — 20 and 
Pr — 0.71. The low-order model predicts a unique stable steady solution for small values of Gr and multiple 
stable steady solutions for large values of Gr in agreement with the direct numerical simulation results based 
on the full model. However, the critical values of Gr at which the multiple stable solutions emerge or lose 
their stability are at variance with the PDE-based calculations. More research is needed in documenting the 
multiplicity of solutions of the PDEs before a meaningful quantitative comparison with the ODE predictions 
can be accomplished. This is a demanding task since direct numerical studies based on the full model are 
computationally intensive and are complicated by the presence of hysteresis effects. 


CONCLUSIONS 

A low- dimensional dynamical model for transitional buoyancy-driven flow in a differentially heated en- 
closure of aspect ratio A=20 has been presented. Empirical eigenfunctions have been determined by ap- 
plying the snapshot version of the proper orthogonal decomposition at “design” conditions Pr Q — 0.71 
and Gr 0 = 3.2 x 10 4 . The computed eigenfunctions are centro-symmetric. Using the four most energetic 
eigenmodes of the velocity field and the four most energetic eigenmodes of the temperature field, a Galerkin 
procedure leads to an eight-equation nonlinear dynamical model. Close to the “design” conditions, the low- 
order model predictions are in excellent agreement with the predictions of the full model. Conditions at the 
onset of the first temporal instability (Hopf bifurcation) of the flow are in excellent agreement with the full 
model predictions in terms of the critical Grashof number, the frequency of oscillations, and the amplitude 
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of oscillations at supercritical conditions. Far from the “design” conditions, the LOM predicts a unique 
steady solution at small values of Gr, and multiple stable steady solutions at large values of Gr, and exhibits 
hysteretic behavior that is qualitatively similar to that observed in direct numerical simulations based on 
the full model. It is believed that low-order models have the potential of becoming a viable tool in the study 
of complex transitional flows. Note, however, that the governing equations have been scaled appropriately 
in order to widen the range of applicability of the resulting low-order model. The possibility of improving 
the performance of the LOM by constructing empirical eigenfunctions that take into account the behavior 
of the flow system at multiple values of Gr is under investigation. 
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Table I 

a) Velocity eigenvalues b) Temperature eigenvalues 


Modes 

Normalized 

Eigenvalue 

Cumulative Energy 
Contribution (%) 

1 

0.6307 

63.07 

2 

0.3163 

94.70 

3 

0.0345 

98.15 

4 

0.0168 

99.83 

5 

0.881xl0 -3 

99.92 

6 

0.677xl0 -3 

99.99 

7 

0.696X10- 4 

99.99 

8 

0.524x 10 -4 

100.0 


Modes 

Normalized 

Eigenvalue 

Cumulative Energy 
Contribution (%) 

i 

0.5583 

55.83 ' 

2 

0.3424 

90.07 

3 

0.0658 

96.65 i 

4 

0.0317 

99.83 

5 

0.921x 10 -3 

99.92 

6 

0.707xl0~ 3 

99.99 

7 

0.603X10" 4 

99.99 

8 

0.293X10' 4 

100.0 
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Grashof number 


Figure 2. Fixed points of the LOM and their stability. Pr = 0.71. 
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